clearvars;
% 常数
global hplank 
global Clight 
global lamdap 
global lamdal 
global Ep 
global El 
global R 
global L 
global delta0
global slope 
global omegap0
global temb 
global z00
global tema 
global H  
global TOC 
global alpha 
global tao1 
global tao2
global sigmae
global gamma1
global gamma2
global nindex 
global fa
global fb
global Leff
global Kc
global dndT
global kesi0
global omegapa
global ft
global omegalaser
global Rp
global PHI
global rp1C
global phi00C
global N0 
global r
global z
global Pout 
global Pin
global g
Pout=1;
hplank=6.626 *10^(-34);% plank常数
Clight=3*10^8;% 光速
lamdap=976*10^(-9);% 泵浦波长
lamdal=2940*10^(-9);% 激光波长
Ep=hplank*Clight/lamdap; % 泵浦光子能量
El=hplank*Clight/lamdal; % 激光光子能量
R =1.5*10^(-3); % 晶体横截面尺寸
L=4*10^(-3);
omega=1.86;% 松弛因子
delta0=0.05;% 内腔损耗 (估计)
omegap0=50*10^(-6); % 聚焦处泵浦光斑半径
slope=0.06; % 泵浦光发散斜率
z00=L/3; % 泵浦聚焦位置
temb=288; % 晶体侧面，热沉温度
tema=298; % 空气温度
H=10; % 空气，晶体热交换系数
Pin=5; % 一端注入的泵浦功率
TOC=0.025;% 输出耦合率
% part 1 不考虑温度影响 （晶体内部为均匀室温），计算N， ，pout, kesi
% 光谱参数
alpha=1400;% 吸收系数
tao1=7250*10^(-6);% 下能级荧光寿命
tao2=120*10^(-6);% 上能级荧光寿命
sigmae =3*10^(-24); % 辐射截面
% 物理参数
fb=0.218;% 上能级玻尔兹曼占有因子
fa=0.055;% 下能级玻尔兹曼占有因子
gamma1=1.3*10^(-21);% 下能级上转换率
gamma2=3.7*10^(-21);% 上能级上转换率
nindex=1.79;% 折射率
Leff=(nindex-1)*L+0.501;% 谐振腔光学长度
Kc=7.6; % 热传导率
dndT=7.3*10^(-6);% 热光系数
kesi0=0.668; % 热沉积比 （不考虑ETU,ESA）
omegapa=(omegap0*L+0.5*slope*((z00)^2+(L-z00)^2)/L);% 晶体内平均泵浦光斑半径
N0=7*10^(15);% 基态粒子数
% STEP1初始温度分布
% 初始热焦距
ft=pi*Kc*omegapa^2/(kesi0*2*Pin*(1-exp (-alpha*L))*dndT);
% 初始激光光斑半径
omegalaser=0.000302916*sqrt(abs((4+49*(1-4/ft)+900*(1+1/100*(-4-49*(1-4/ft))-4/ft))/sqrt(4-(51/100*(1+900*(1/100*(-1+49/ft)-1/ft)-49/ft)+(1-49/ft)*(1+1/100*(-4-49*(1-4/ft))-4/ft))^2))); % 晶体处激光光斑半径-单位米
m=(omegalaser/omegapa)^2;
Rp=Pin*(1-exp (-alpha*L))/Ep; % 泵浦速率
% 求解场分布函数中的常数部分
rp1C=2*alpha/(pi*omegapa^2* (1-exp (-alpha*L))); % 泵浦场的常数部分
phi00C=2*nindex/(pi*omegalaser^2*Leff);% 激光场的常数部分
% 激光器输出功率
PHI=2*Leff*Pout/((-log(1-TOC)*El*Clight));
B=nindex*(delta0-log(1-TOC))/(2*Leff*sigmae);%等式右端
N1=(-1+sqrt (1+4*gamma1*Rp*N0*tao1^2))/2*gamma1*tao1;% 可以先积分出来
B1=(sigmae*Clight*fb^2/(gamma2*nindex*pi*Leff*omegalaser^2))^2;% phi平方前的常数
B2=(sigmae*Clight*fb^3/(tao2*gamma2*nindex*pi*Leff*omegalaser^2));% phi前的常数
B3=(fb/2*tao2*gamma2)^2;
B4=(sigmae*Clight*fa*fb^2/(tao1*gamma1*gamma2*nindex*pi*Leff*omegalaser^2));% phi前的常数
B5= fb^2/(2*gamma2*gamma1*tao1^2);
B6=sqrt (4*gamma1*tao1^2*N0*Pin*2*alpha/(Ep*pi*omegapa^2));
B7=2*fb^2*N0* Pin*2*alpha/ (gamma2*Ep*pi*omegapa^2);% 以上为等式左端第一项的常数
B8=fb /(tao2*gamma2*pi* omegalaser^2*Leff);
B9=2*sigmae *fb^2/( nindex*gamma2* Leff *pi^2*omegalaser^4*(-log(1-TOC))*El);
B10=2*pi*L*(-omegalaser^2/4);
a=(B8+B9*Pout)*B10;% 以上为左端第二项 (已经是最终积分结果)
B11=fa/(tao1*gamma1* pi* omegalaser^2*Leff)*(-omegalaser^2/4);% 常数
B12=sqrt (8* tao1^2*gamma1*alpha*N0/(Ep* pi* omegapa^2));% pin前边为常数
B13=B11*2*pi*L;% 常数
F3=sqrt (exp (-alpha*z)*exp (-2*r^2/omegapa^2)+1)*exp (-2*r^2/omegalaser^2);% 积分=0.0042
F4=B11*B12* F3*sqrt(Pin)-B13;% 含有pin的关系式,pin在第一项
% 左端第三项 （只需将F3积分）
 f=2/(pi* omegalaser.^2*Leff)*sqrt(B1*PHI.^2*exp(-8*r.^2/omegalaser.^2)+B2* PHI* exp(-6*r.^2/omegalaser.^2)+B3*exp(-4*r.^2/omegalaser.^2)+exp(-4*r.^2/omegalaser.^2)*(B5-B4* PHI* exp(-2*r.^2/omegalaser.^2))*(1-sqrt(B6* exp (-alpha*z)*exp (-2*r.^2/omegapa.^2)+1))+B7* exp (-alpha*z)*exp (-2*r.^2/omegapa.^2));
